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Abstract 

We determine the topological susceptibility of the gauge configurations generated by lattice 
simulations using two flavors of optimal domain-wall fermion on the 16'^ x 32 lattice with length 
16 in the fifth dimension, at the lattice spacing a ~ 0.1 fm. Using the adaptive thick-restart 
Lanczos algorithm, we project the low-lying eigenmodes of the overlap Dirac operator, and obtain 
the topological charge of each configuration, for eight ensembles with pion masses in the range 
220 — 550 MeV. From the topological charge, we compute the topological susceptibility and the 
second normalized cumulant. Our result of the topological susceptibility agrees with the sea- 
quark mass dependence predicted by the chiral perturbation theory and provides a determination 
of the chiral condensate, ^'^^(2 GeV) = [259(6) (7) MeYf, and the pion decay constant F-^ = 
92(12)(2) MeV. 
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I. INTRODUCTION 



The vacuum of Quantum Chromodynamics (QCD) has a non-trivial topological struc- 
ture. The cluster property and the gauge invariance require that the ground state must be 
the 6 vacuum, a superposition of gauge configurations in different topological sectors. The 
topological susceptibility (xt) is the most crucial quantity to measure the topological fluctu- 
ations of the QCD vacuum, which plays an important role in breaking the Ua{^) symmetry. 
Theoretically, Xt is defined as 

Xt = j (fx (p(a;)p(0)) , (1) 

where p{x) = e^,yAcrtr[F^,^(x)FAo-(x)]/(327r^), is the topological charge density expressed in 
term of the matrix- valued field tensor F^^. With mild assumptions, Witten [1] and Veneziano 
[2I obtained a relationship between the topological susceptibility in the quenched approxi- 
mation and the mass of rj' meson (flavor singlet) in the full QCD. This implies that the mass 
of rj' is essentially due to the axial anomaly relating to non-trivial topological fluctuations, 
unlike those of the (non-singlet) approximate Goldstone bosons. 
From ([1]), we obtain 



Xt 



Qt = J d^xp{x), (2) 



where Q is the volume of the system, and Qt is the topological charge (which is an integer 
for QCD). Thus, one can determine Xt by counting the number of gauge configurations for 
each topological sector. Furthermore, we can also obtain the second normalized cumulant 

1 



{Qt)-my , (3) 



which is related to the leading anomalous contribution to the rj' — rj' scattering amplitude 
in QCD, as well as the dependence of the vacuum energy on the vacuum angle 6. (For a 
recent review, see, for example, Ref. 3|] and references therein.) 

However, for lattice QCD, it is difficult to extract p{x) and Qt unambiguously from the 
gauge link variables, due to their rather strong fluctuations. To circumvent this difficulty, 
we may consider the Atiyah-Singer index theorem Qt = — n_ = index(P), where 
n± is the number of zero modes of the massless Dirac operator V = 7/^(9^ + igA^j^) with ± 
chirality. 
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For lattice QCD with exact chiral symmetry, it is well-known that the overlap Dirac 
operator js, 6| in a topologically non-trivial gauge background possesses exact zero modes 
(with definite chirality) satisfying the Atiyah-Singer index theorem. Thus we can obtain the 
topological charge from the index of the overlap Dirac operator. Writing the overlap Dirac 
operator as 



ITlQ 1+75- 



where is the standard Hermitian Wilson operator with negative mass —rriQ (0 < mo < 2), 
then its index is 



index (D) = Tr 



75 1 



Do \ 
2moJ 



where Tr denotes trace over Dirac, color, and site indices. 

In this paper, we measure the topological charge of the gauge configurations generated by 
lattice simulations of two fiavors QCD on a 16^ x 32 lattice, with the optimal domain-wall 
fermion (ODWF) [3| at A^^ = 16, and plaquette gauge action at /3 = 5.95, for eight sea-quark 
masses m^a = 0.01, ■ ■ ■ , 0.08 with the interval 0.01. 

Mathematically, ODWF is a theoretical framework which can preserve the chiral sym- 
metry optimally for any given Ng, with a set of analytical weights {ws, s = 1, ■ ■ ■ , N^}, one 
for each layer in the fifth dimension I?]. Thus the artifacts due to the chiral symmetry 
breaking with finite Ng can be reduced to the minimum, especially in the chiral regime. The 
4-dimensional effective Dirac operator of massless ODWF is 



i-nf=iT, 



D = mo[l + 755'opt(-ff«;)], Sopt{Hu 



Tg 



ojgH,. 



(4) 



which is exactly equal to the Zolotarev optimal rational approximation of the overlap Dirac 
operator. That is, SrmtiHw) = H^RziHw), where Rz{Hw) is the optimal rational approxi- 



mation of {Hi 



-1/2 



a. 



We use the adaptive thick-restart Lanczos algorithm 10| to project the low-lying eigen- 
modes of the 4-dimensional effective Dirac operator and obtain the topological charge 
Qt of each gauge configuration. Then we compute the topological susceptibilty Xt and the 
second normalized cumulant C4, and compare our results to the Chiral Perturbation Theory 
(ChPT). We summarize the ChPT formulas as follows. 
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In 1992, Leutwyler and Smilga ll| derived the relationship between xt and the quark 



mass, at the leading order in ChPT. For 2 flavors QCD, it reads 

Xt = S (m-i + m^^) ^ , (5) 

where m^, are the quark masses, and E is the chiral condensate. This implies that in 
the chiral limit (m„ — t- 0) the topological susceptibility is suppressed by the internal quark 
loops. Most importantly, (|5]) provides a viable way to extract S from Xt in the chiral regime. 

Recently, the topological suscept^ility has been derived to the one-loop order in ChPT 
for an arbitrary number of flavors [12]. For Nj = 2 with degenerate u and d quark masses 
{niu = rrid = nig), the foumula reduces to 



niq 2 



where Li are renormalized low-energy coupling constants deflned at figub [l3|. In this paper, 
we flx fisub = 770 MeV. 

Furthermore, the second normalized cumulant Ca has been derived in ChPT at the tree- 
level for an arbitrary number of flavors . For 2 flavors QCD, it reads 

C4 = -S {mZ' + m/) {m-' + m, ^) . (7) 

In the isospin limit {rriu = m^), the ratio Ci/xt goes to —1/4. 

In this paper, we investigate to what extent the quark mass dependences of Xt in lattice 
QCD with optimal domain-wall fermion would agree with the ChPT to the one-loop order, 
and to determine S and from our data of Xt- In principle, we can also extract S from C4, 
however, this would require much higher statistics than that of Xt- 



II. LATTICE SETUP 



Simulations are carried out for two flavors QCD on a 16^ x 32 lattice at the lattice spacing 
a ~ 0.1 fm, for eight sea-quark masses m^a = 0.01,0.02,0.03,0.04,0.05,0.06,0.07, and 0.08 
respectively. For the quark part, we use the optimal domain-wall fermion with A^^ = 16. 
For the gluon part, we use the plaquette action at /3 = 5.90 and (3 = 5.95 respectively. 
An outline of our simulation algorithm and its acceleration with Nvidia CPUs has been 
p.sen.d . Re.. Q Q. a.d the detai. wUl be p...ed in Ref. Q. Ou. p.elin,ina. 



physical results of the /3 = 5.90 ensemble have been presented in Refs. jl7, ll8|]. In this 
paper, we present our results of the topological susceptibility Xt and the second normahzed 
cumulant C4 of the (3 = 5.95 ensemble. 

For each sea-quark mass, we perform hybrid Monte-Carlo simulations on 30 CPUs inde- 
pendently, with each GPU generating 400 trajectories. After discarding 300 trajectories for 
thermalization, each GPU yields 100 trajectories. Thus, with 30 GPUs running indepen- 
dently, we accumulated total 3000 trajectories for each sea-quark mass. From the saturation 
of the binning error of the plaquette, as well as the evolution of the topological charge, we 
estimate the autocorrelation time to be around 10 trajectories. Thus we sample one con- 
figuration every 10 trajectories, then we have 300 configurations for each sea-quark mass. 
With a GPU cluster of 250 GPUs, we can simulate 8 sea-quark masses concurrently. It takes 
about 5 months to complete the simulations for the /3 = 5.95 ensemble. 

We determine the lattice spacing by heavy quark potential with Sommer parameter 
ro = 0.49 fm. Using the linear fit, we obtain the lattice spacing in the chiral limit, 
a = 0.1032(2) fm, which gives a^^ = 1.911(4)(6) GeV, where the systematic error is es- 
timated with the uncertainty of tq. 

For each configuration, we calculate the zero modes plus 80 conjugate pairs of the lowest- 
lying eignmodes of the overlap Dirac operator. We outline our procedures as follows. First, 
we project 240 low-lying eigenmodes of using adaptive thick-restart Lanczos alogorithm 
(a-TRLan) lOj, where each eigenmode has a residual less than 10~^^. Then we approximate 
the sign function of the overlap operator by the Zolotarev optimal rational approximation 
with 64 poles, where the coefficents are fixed with X'^^x = (6-4)^, and A^j„ equal to the 
maximum of the 240 projected eigenvalues of H^. Then the sign function error is less than 
10~^^. Using the 240 low-modes of and the Zolotarev approximation with 64 poles, we 
use the a-TRLan algorithm again to project the zero modes plus 80 conjugate pairs of the 
lowest-lying eignmodes of the overlap operator, where each eigenmode has a residual less 
than 10~^^. We store all projected eigenmodes for the later use. In this paper, we use the 
index of the zero modes to compute Xt and C4. 
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FIG. 1: Histogram of topological charge distribution for eight sea quark masses, m^a 
0.01,0.02,0.03,0.04,0.05,0.06,0.07, and 0.08 respectively. 



III. RESULTS 

In Fig. [H we plot the histogram of topological charge distribution for m^a = 
0.01, 0.02, ■ ■ ■ , 0.08 respectively. Evidently, the probability distribution of Qt for each sea- 
quark mass behaves like a Gaussian, and it becomes more sharply peaked around Qt = as 
the sea-quark mass nig gets smaller. 

Using the result of Qt, we compute the topological susceptibility Xt and the second 
normalized cumulant C4 ([3]). In Table HI we list our results of Xt, ca, and the ratio c^/xt- 
The error is estimated using the jackknife method with bin size of 10 configurations, with 
which the statistical error saturates. 

Evidently, the statistical error of the topological susceptibility is about 10%, while that 
of C4 is very large due to low statistics. Therefore, we cannot draw any conclusions from our 
result of C4, as well as from the ratio Ci/xt- 

In Fig [21 we plot our data of Xt versus the sea quark mass m^. The data points of Xt 



6 



niqa 


Xt 




C4 


c^/Xt 


0.01 


1.13(10) X 10" 


-5 


-1.22(1.15) X 10-^ 


-1.07(1.01) 


0.02 


2.24(18) X 10" 


-5 


-3.79(2.77) X 10^5 


-1.69(1.24) 


0.03 


3.29(27) X 10" 


-5 


-1.25(2.87) X 10"^ 


-0.04(0.87) 


0.04 


4.40(30) X 10" 


-5 


6.39(4.16) X 10"^ 


1.45(95) 


0.05 


5.31(41) X 10- 


-5 


5.75(7.96) X 10-5 


1.08(1.50) 


0.06 


6.04(44) X 10" 


-5 


1.00(1.09) X 10-"^ 


1.66(1.81) 


0.07 


7.24(55) X 10- 


-5 


-5.12(134) X 10-6 


-0.07(1.86) 


0.08 


7.01(79) X 10" 


-5 


-7.10(6.58) X 10-^ 


-10.13(9.46) 



TABLE I: The topological susceptibility xt, the second normalized cumulant C4, and their ratio 
Ci/xt, versus the sea quark masses, for Nj = 2 lattice QCD with the optimal domain-wall fermion. 
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FIG. 2: The topological susceptibility xt versus the sea-quark mass niq for 2 flavors lattice QCD 
with ODWF. The straight line is the fit with the LO ChPT ([5]). 

are well fitted by the Leutwyler and Smilga formula ([5]) with Sa^ = 0.00200(15). The fitted 
curve is plotted as the solid line in Fig El 

In Fig. [3l we plot our data of Xt/^q versus the sea-quark mass rriq. The data points of 
Xt/rrig are well fitted by the NLO ChPT formula dSD with Sa^ = 0.0020(2), F^a = 0.048(7), 
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FIG. 3: The ratio Xt/'mg versus the sea-quark mass mq for 2 flavors lattice QCD with ODWF. The 
solid line is the fit with the NLO ChPT dH). 



and 



(2L6 + 2L7 + L8) = -0.0001(3), 



(8) 



where fisub = 770 MeV has been used. Using a ^ = 1.911(4)(6) GeV, we obtain S 
[241(6) (1) MeV]3, and 

= 92(12) (2) MeV, 



(9) 



where the errors represent a combined statistical error and the systematic error respectively. 

In order to convert S to that in the MS scheme, we calculate the renormalization factor 
Zf^^{2 GeV) using the non-perturbative renormalization technique through the RI/MOM 



scheme [2^, and our result is [21] 

Zf^{2 GeV) = 1.244(18)(39). 
Then the value of E is transcribed to 

E^(2 GeV) = [259(6) (7) MeV]^ 



(10) 



where the errors represent a combined statistical error (a~^ and Zf^^) and the systematic 
error respectively. Since the present calculation is done at a single lattice spacing, the 
discretization error cannot be quantified reliably, but we do not expect much larger error 
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because the optimal domain- wall fermion action is free from 0(a) discretization effects. Our 
result of E (|T0|) is in good agreement with that extracted from xt in (2+1) flavors QCD with 



domain-wall fermion 



22|, as well as with those extracted from Xt in ^/ = 2 and Nf = 2 + 1 



lattice QCD with overlap fermion in a fixed topology 



23 



24|. 



IV. CONCLUDING REMARK 



To summarize, we measure the topological charge of the gauge configurations generated 
by lattice simulations of 2 fiavors QCD with the optimal domain-wall fermion at A^^ = 16 and 
plaquette gauge action at /3 = 5.95, on a 16'^ x 32 lattice. We use the adaptive thick-restart 
Lanczos algorithm to compute the low-lying eigenmodes of the overlap Dirac operator, and 
obtain the topological charge of each gauge configuration, and from which we compute the 
topological susceptibility for 8 sea-quark masses, each of 300 configurations. Our result of 
the topological susceptibility agrees with the sea-quark mass dependence predicted by the 
NLO ChPT formula ([H]), and gives the first determination of both the pion decay constant ^ 
and the chiral condensate (fTOj) simultaneously from the topological susceptibility. Moreover, 
this study shows that it is feasible to perform large-scale simulations of unquenched lattice 
QCD, which not only preserve the chiral symmetry to a high precision, but also sample all 
topological sectors ergodically. 
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